/*
 * vectorvolume.cpp
 *
 * Copyright 2012 The University of Auckland
 *
 *  Created on: 28/08/2012
 *      Author: cvan680
 */

#include "vectorvolume.h"
#include "voreen/core/datastructures/volume/volumefactory.h"
#include "voreen/core/datastructures/volume/volumeatomic.h"
#include "voreen/core/datastructures/volume/volumeram.h"
#include "voreen/core/datastructures/volume/volume.h"
#include "tgt/textureunit.h"
#include "tgt/filesystem.h"

#include <Eigen/Core>
#include <Eigen/Eigen>

namespace voreen {

VectorVolume::VectorVolume()
	: VolumeProcessor()
	, outport(Port::OUTPORT, "volume", 0)
	, prop_xdim("dimension_x", "Dimension (x)", 128, 1, 1024)
	, prop_ydim("dimension_y", "Dimension (y)", 128, 1, 1024)
	, prop_zdim("dimension_z", "Dimension (z)", 2, 1, 1024)
, prop_sing1Type(0)
, prop_sing1_x("sing1_x", "Singularity 1 location x", 50, 0, 100)
, prop_sing1_y("sing1_y", "Singularity 1 location y", 50, 0, 100)
, prop_sing2Type(0)
, prop_sing2_x("sing2_x", "Singularity 2 location x", 50, 0, 100)
, prop_sing2_y("sing2_y", "Singularity 2 location y", 50, 0, 100)
, prop_sing3Type(0)
, prop_sing3_x("sing3_x", "Singularity 3 location x", 50, 0, 100)
, prop_sing3_y("sing3_y", "Singularity 3 location y", 50, 0, 100)
, prop_sing4Type(0)
, prop_sing4_x("sing4_x", "Singularity 4 location x", 50, 0, 100)
, prop_sing4_y("sing4_y", "Singularity 4 location y", 50, 0, 100)
, prop_sing5Type(0)
, prop_sing5_x("sing5_x", "Singularity 5 location x", 50, 0, 100)
, prop_sing5_y("sing5_y", "Singularity 5 location y", 50, 0, 100)
, prop_sing6Type(0)
, prop_sing6_x("sing6_x", "Singularity 6 location x", 50, 0, 100)
, prop_sing6_y("sing6_y", "Singularity 6 location y", 50, 0, 100)
	, prop_decay("decay","Decay", 0.001f, 0.0f, 1.0f, Processor::VALID)
	, prop_useFirstEigenvector("useFirstEigenVector", "Use first eigenvector", true)
	, update_prop("update", "Update outputs")
{
    addPort(outport);

    addProperty(prop_xdim);
    addProperty(prop_ydim);
    addProperty(prop_zdim);

    prop_sing1Type = new OptionProperty<SingularityType>("singularityType1", "Singularity 1: ");
    prop_sing1Type->addOption("none", "----", NONE);
    prop_sing1Type->addOption("node", "node", NODE);
    prop_sing1Type->addOption("saddle", "saddle", SADDLE);
    prop_sing1Type->addOption("center", "center", CENTER);
    prop_sing1Type->addOption("attr_focus", "attr_focus", ATTRACTING_FOCUS);
    prop_sing1Type->addOption("rep_focus", "rep_focus", REPELLING_FOCUS);
	addProperty(prop_sing1Type);
	addProperty(prop_sing1_x);
	addProperty(prop_sing1_y);

    prop_sing2Type = new OptionProperty<SingularityType>("singularityType2", "Singularity 2: ");
    prop_sing2Type->addOption("none", "----", NONE);
    prop_sing2Type->addOption("node", "node", NODE);
    prop_sing2Type->addOption("saddle", "saddle", SADDLE);
    prop_sing2Type->addOption("center", "center", CENTER);
    prop_sing2Type->addOption("attr_focus", "attr_focus", ATTRACTING_FOCUS);
    prop_sing2Type->addOption("rep_focus", "rep_focus", REPELLING_FOCUS);
	addProperty(prop_sing2Type);
	addProperty(prop_sing2_x);
	addProperty(prop_sing2_y);

    prop_sing3Type = new OptionProperty<SingularityType>("singularityType3", "Singularity 3: ");
    prop_sing3Type->addOption("none", "----", NONE);
    prop_sing3Type->addOption("node", "node", NODE);
    prop_sing3Type->addOption("saddle", "saddle", SADDLE);
    prop_sing3Type->addOption("center", "center", CENTER);
    prop_sing3Type->addOption("attr_focus", "attr_focus", ATTRACTING_FOCUS);
    prop_sing3Type->addOption("rep_focus", "rep_focus", REPELLING_FOCUS);
	addProperty(prop_sing3Type);
	addProperty(prop_sing3_x);
	addProperty(prop_sing3_y);

    prop_sing4Type = new OptionProperty<SingularityType>("singularityType4", "Singularity 4: ");
    prop_sing4Type->addOption("none", "----", NONE);
    prop_sing4Type->addOption("node", "node", NODE);
    prop_sing4Type->addOption("saddle", "saddle", SADDLE);
    prop_sing4Type->addOption("center", "center", CENTER);
    prop_sing4Type->addOption("attr_focus", "attr_focus", ATTRACTING_FOCUS);
    prop_sing4Type->addOption("rep_focus", "rep_focus", REPELLING_FOCUS);
	addProperty(prop_sing4Type);
	addProperty(prop_sing4_x);
	addProperty(prop_sing4_y);

    prop_sing5Type = new OptionProperty<SingularityType>("singularityType5", "Singularity 5: ");
    prop_sing5Type->addOption("none", "----", NONE);
    prop_sing5Type->addOption("node", "node", NODE);
    prop_sing5Type->addOption("saddle", "saddle", SADDLE);
    prop_sing5Type->addOption("center", "center", CENTER);
    prop_sing5Type->addOption("attr_focus", "attr_focus", ATTRACTING_FOCUS);
    prop_sing5Type->addOption("rep_focus", "rep_focus", REPELLING_FOCUS);
	addProperty(prop_sing5Type);
	addProperty(prop_sing5_x);
	addProperty(prop_sing5_y);

    prop_sing6Type = new OptionProperty<SingularityType>("singularityType6", "Singularity 6: ");
    prop_sing6Type->addOption("none", "----", NONE);
    prop_sing6Type->addOption("node", "node", NODE);
    prop_sing6Type->addOption("saddle", "saddle", SADDLE);
    prop_sing6Type->addOption("center", "center", CENTER);
    prop_sing6Type->addOption("attr_focus", "attr_focus", ATTRACTING_FOCUS);
    prop_sing6Type->addOption("rep_focus", "rep_focus", REPELLING_FOCUS);
	addProperty(prop_sing6Type);
	addProperty(prop_sing6_x);
	addProperty(prop_sing6_y);


	addProperty(prop_decay);

	addProperty(prop_useFirstEigenvector);

	update_prop.onClick(CallMemberAction<VectorVolume>(this, &VectorVolume::update));
	addProperty(update_prop);
}


Processor* VectorVolume::create() const {
    return new VectorVolume();
}

/**
 * What to do here? Now we do everything through the update button.
 * We don't want to create new volumes everytime we change any parameter,
 * so just leave this empty and only update when user presses button.
 */
void VectorVolume::process() {

    // check for OpenGL errors
    LGL_ERROR;
}


void VectorVolume::initialize() throw (VoreenException) {
    // call superclass function first
	VolumeProcessor::initialize();
	volumefactory = new VolumeFactory();
}

void VectorVolume::deinitialize() throw (VoreenException) {
    // call superclass function last
	VolumeProcessor::deinitialize();
}


void VectorVolume::update() {
	volume_ = multipleSingularitiesVolume();
	outport.setData(volume_, true);
}


Volume* VectorVolume::saddleVectorVolume(){

	tgt::ivec3 dims = getDimensions();

	// 1. Calculate Tensor field with saddle
	float v11 = 1.0;
	float v12 = 0.0;
	float v21 = 0.0;
	float v22 = 1.0;

	tgt::ivec2 saddle(prop_sing1_x.get(), prop_sing1_y.get());
	float dist = 0.0;

	//TODO is this meant to be between 0 and 1?
	// decay	factor at dist=100
	// .00010	.368
	// .00025	.082
	// .00050	.007
	float decay = prop_decay.get();
	float factor = 0.0;

	// 2. calculate eigenvectors per tensor
	VolumeRAM* vectorVR = volumefactory->create("Vector3(float)", tgt::svec3(dims.x, dims.y, dims.z));
	tgt::vec2 vect(0.0, 0.0);
	std::cout << "\n\n";
	for(int x=0; x<dims.x; x++){
		for(int y=0; y<dims.y; y++){
			//vect = tfield->get(x,y).eigenVector3D();
			tgt::ivec2 v(saddle.x-x, saddle.y-y);
			dist = sqrtf(tgt::dot(v, v));
			factor = exp(-1 * decay * (dist*dist));

			switch(prop_sing1Type->getValue()){
				case ATTRACTING_FOCUS: //TODO these two are not configured
					// [ a	b ]
					// [ c	d ]
					v11 = factor * x;			// a
					v12 = factor * -1.0 * y;	// b
					v21 = factor * -1.0 * y;	// c
					v22 = factor * -1.0 * x;	// d
					break;
				case REPELLING_FOCUS:
					// [ a	b ]
					// [ c	d ]
					v11 = factor * x;			// a
					v12 = factor * -1.0 * y;	// b
					v21 = factor * -1.0 * y;	// c
					v22 = factor * -1.0 * x;	// d
					break;
				case NODE:
					v11 = factor * (x*x - y*y);
					v12 = factor * (2.0*x*y);
					v21 = factor * (2.0*x*y);
					v22 = factor * -1.0 * (x*x - y*y);
					break;
				case SADDLE:
					v11 = factor * (x*x - y*y);
					v12 = factor * (-2.0*x*y);
					v21 = factor * (-2.0*x*y);
					v22 = factor * -1.0 * (x*x - y*y);
					break;
				case CENTER:
					v11 = factor * (y*y - x*x);
					v12 = factor * (-2.0*x*y);
					v21 = factor * (-2.0*x*y);
					v22 = factor * -1.0 * (y*y - x*x);
					break;
			}
			bool show = false;// x==y && x < int(dims.x*0.55) && x > int(dims.x*0.45);

			vect = calcEigenvector(tgt::vec2(v11, v21), tgt::vec2(v12, v22), show);
			for(int z=0; z<dims.z; z++){
				vectorVR->setVoxelNormalized(vect.x, x, y, z, 0);
				vectorVR->setVoxelNormalized(vect.y, x, y, z, 1);
				vectorVR->setVoxelNormalized(0.0, 	 x, y, z, 2);
			}
		}
	}

	Volume* result = new Volume(vectorVR, tgt::vec3(1.0, 1.0, 1.0), tgt::vec3(0.0, 0.0, 0.0));
	centerVolume(result);
	return result;
}



Volume* VectorVolume::synthesizeVectorVolume(){

	tgt::ivec3 dims = getDimensions();

	tgt::ivec2 singularity(prop_sing1_x.get(), prop_sing1_y.get());
	float dist = 0.0;

	//TODO is this meant to be between 0 and 1?
	// decay	factor at dist=100
	// .00010	.368
	// .00025	.082
	// .00050	.007
	float decay = prop_decay.get();
	float factor = 0.0;

	float xvec = 0.0;
	float yvec = 0.0;

	// 2. calculate eigenvectors per tensor
	VolumeRAM* vectorVR = volumefactory->create("Vector3(float)", tgt::svec3(dims.x, dims.y, dims.z));
	tgt::vec2 vect(0.0, 0.0);
	for(int x=0; x<dims.x; x++){
		for(int y=0; y<dims.y; y++){
			//vect = tfield->get(x,y).eigenVector3D();
			tgt::ivec2 v(singularity.x-x, singularity.y-y);
			dist = sqrtf(tgt::dot(v, v));
			factor = exp(-1 * decay * (dist*dist));

			switch(prop_sing1Type->getValue()){
				case ATTRACTING_FOCUS: //TODO these two are not configured
					xvec = -v.y - .5*v.x;
					yvec = v.x - .5*v.y;
					break;
				case REPELLING_FOCUS: //TODO these two are not configured
					xvec = -v.y + .5*v.x;
					yvec = v.x + .5*v.y;
					break;
				case NODE:
					xvec = v.x+3*v.y;
					yvec = -2*v.x-4*v.y;
					break;
				case SADDLE:
					xvec = v.x;
					yvec = -v.y;
					break;
				case CENTER:
					xvec = v.x+-4*v.y;
					yvec = -4*v.x + v.y;
					break;
			}

			for(int z=0; z<dims.z; z++){
				vectorVR->setVoxelNormalized(xvec, x, y, z, 0);
				vectorVR->setVoxelNormalized(yvec, x, y, z, 1);
				vectorVR->setVoxelNormalized(0.0,  x, y, z, 2);
			}
		}
	}

	Volume* result = new Volume(vectorVR, tgt::vec3(1.0, 1.0, 1.0), tgt::vec3(0.0, 0.0, 0.0));
	centerVolume(result);
	return result;
}


Volume* VectorVolume::multipleSingularitiesVolume(){
	tgt::ivec3 dims = getDimensions();

	tgt::vec2 singularity1(dims.x * (prop_sing1_x.get()/100.0), dims.y * prop_sing1_y.get()/100.0);
	tgt::vec2 singularity2(dims.x * (prop_sing2_x.get()/100.0), dims.y * prop_sing2_y.get()/100.0);
	tgt::vec2 singularity3(dims.x * (prop_sing3_x.get()/100.0), dims.y * prop_sing3_y.get()/100.0);
	tgt::vec2 singularity4(dims.x * (prop_sing4_x.get()/100.0), dims.y * prop_sing4_y.get()/100.0);
	tgt::vec2 singularity5(dims.x * (prop_sing5_x.get()/100.0), dims.y * prop_sing5_y.get()/100.0);
	tgt::vec2 singularity6(dims.x * (prop_sing6_x.get()/100.0), dims.y * prop_sing6_y.get()/100.0);

	int s1type = prop_sing1Type->getValue();
	int s2type = prop_sing2Type->getValue();
	int s3type = prop_sing3Type->getValue();
	int s4type = prop_sing4Type->getValue();
	int s5type = prop_sing5Type->getValue();
	int s6type = prop_sing6Type->getValue();

	tgt::vec2 value(0.0, 0.0);
	tgt::vec2 v1(0, 0);
	tgt::vec2 v2(0, 0);
	tgt::vec2 v3(0, 0);
	tgt::vec2 v4(0, 0);
	tgt::vec2 v5(0, 0);
	tgt::vec2 v6(0, 0);

	// 2. calculate eigenvectors per tensor
	VolumeRAM* vectorVR = volumefactory->create("Vector3(float)", tgt::svec3(dims.x, dims.y, dims.z));
	tgt::vec2 vect(0.0, 0.0);
	for(int x=0; x<dims.x; x++){
		for(int y=0; y<dims.y; y++){
			value = tgt::vec2(0.0, 0.0);

			v1 = tgt::vec2(singularity1.x-x, singularity1.y-y);
			v2 = tgt::vec2(singularity2.x-x, singularity2.y-y);
			v3 = tgt::vec2(singularity3.x-x, singularity3.y-y);
			v4 = tgt::ivec2(x-singularity4.x, y-singularity4.y);
			v5 = tgt::ivec2(x-singularity5.x, y-singularity5.y);
			v6 = tgt::ivec2(x-singularity6.x, y-singularity6.y);

			if(s1type != NONE){ value += calcVec(v1, s1type); }
			if(s2type != NONE){ value += calcVec(v2, s2type); }
			if(s3type != NONE){ value += calcVec(v3, s3type); }
			if(s4type != NONE){ value += calcVec(v4, s4type); }
			if(s5type != NONE){ value += calcVec(v5, s5type); }
			if(s6type != NONE){ value += calcVec(v6, s6type); }

			for(int z=0; z<dims.z; z++){
				vectorVR->setVoxelNormalized(value.x, x, y, z, 0);
				vectorVR->setVoxelNormalized(value.y, x, y, z, 1);
				vectorVR->setVoxelNormalized(0.0,  	  x, y, z, 2);

			}
		}
	}

	Volume* result = new Volume(vectorVR, tgt::vec3(1.0, 1.0, 1.0), tgt::vec3(0.0, 0.0, 0.0));
	centerVolume(result);
	return result;
}


tgt::vec2 VectorVolume::calcVec(tgt::vec2 relcoords, int type){
	float decay = prop_decay.get() * prop_decay.get() * prop_decay.get();
	float dist = sqrtf(tgt::dot(relcoords, relcoords));
	float factor = exp(-1 * decay * (dist*dist));

	if(prop_useFirstEigenvector.get()){
		factor = 1.0; // just in case 'factor' messes things up
	}

	float angular_velocity = 0.15;

	float xval = 0.0;
	float yval = 0.0;

	switch(type){
		case NODE:
			xval = factor * relcoords.x;
			yval = factor * relcoords.y;
			break;
		case SADDLE:
			xval = factor * relcoords.x;
			yval = factor * - relcoords.y;
			break;
		case ATTRACTING_FOCUS:
			xval = factor * angular_velocity * (-1*relcoords.y - .5*relcoords.x);
			yval = factor * angular_velocity * (relcoords.x - .5*relcoords.y);
			break;
		case REPELLING_FOCUS:
			xval = factor * angular_velocity * (-1*relcoords.y + .5*relcoords.x);
			yval = factor * angular_velocity * (relcoords.x + .5*relcoords.y);
			break;
		case CENTER:
			xval = factor * (-1.0 * angular_velocity * relcoords.y);
			yval = factor * (angular_velocity * relcoords.x);
			break;
		case NONE:
		default:
			break;
	}

	return tgt::vec2( xval, yval );
}



float VectorVolume::calculateVectorMagnitude( tgt::vec3 v ) {
	return sqrtf(tgt::dot(v, v));
}

tgt::vec3 VectorVolume::normalizeVector( tgt::vec3 v ) {
	return v / calculateVectorMagnitude(v);
}

/**
 * Helper function to calculate eigenvector of 2x2 matrix with
 * c1 and c2 as first and second column respectively
 */
tgt::vec2 VectorVolume::calcEigenvector(tgt::vec2 c1, tgt::vec2 c2, bool show){

	using namespace Eigen;

	MatrixXf tmp = MatrixXf::Zero(2, 2);

	tmp(0,0) = c1.x;
	tmp(1,0) = c1.y;
	tmp(0,1) = c2.x;
	tmp(1,1) = c2.y;

    EigenSolver<MatrixXf> solver(tmp);

    MatrixXf eigenVectors = MatrixXf::Zero(2, 2);
    eigenVectors = solver.eigenvectors().real();

    //TODO perhaps NOT normalize?

    /*
     * [ a	b ] -->	[ c1.x	c2.x ]
     * [ c	d ]		[ c1.y	c2.y ]
     *
     * eigenvectors: [a c] [b d]
     *
     * (0,0) = a
     * (1,0) = c
     * (0,1) = b
     * (1,1) = d
     */

    tgt::vec2 ev(0.0, 0.0);
    if(prop_useFirstEigenvector.get()){
    	ev = tgt::vec2(eigenVectors(0, 0), eigenVectors(1, 0));
    }else{
    	ev = tgt::vec2(eigenVectors(0, 1), eigenVectors(1, 1));
    }
    return ev;
}

tgt::ivec3 VectorVolume::getDimensions(){
	return tgt::ivec3(prop_xdim.get(), prop_ydim.get(), prop_zdim.get());
}

} // namespace
